##########################################################################################

library(data.table)
library(optparse)
library(ArchR)

##########################################################################################
option_list <- list(
    make_option(c("--motif_all_file"), type = "character"),
    make_option(c("--comine_data_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 整合atac和rna的文件
    motif_all_file <- "~/20231121_singleMuti/results/celltype_plot/peak2gene/germ_mfuzz/Motif.Peak-Gene.positive.rds"

    ## 整合atac和rna的文件
    comine_data_file <- "/public/home/xxf2019/20231121_singleMuti/results/subcell/cluster2/cluster2.combineRNA.motif_peak2gene.Rdata"

    ## 输出
    out_path <- "/public/home/xxf2019/20231121_singleMuti/results/celltype_plot/sperm_recluster/"

    #comine_data_all_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ_peak-gene/testis_combined_peak.combineRNA.qc.Rdata"


}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

comine_data_file <- opt$comine_data_file
motif_all_file <- opt$motif_all_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################

a <- load(comine_data_file)
b <- readRDS(motif_all_file)
#c <- load(comine_data_all_file)

###########################################################################################
proj <- atac_proj

set.seed(123)
## 重聚类
# resolution越大，图上越紧密
proj <- addIterativeLSI(ArchRProj = proj, useMatrix = "TileMatrix", 
    clusterParams = list(resolution = c(0.6), sampleCells = 500, maxClusters = 2, n.start = 10) ,
    name = "IterativeLSI" , force = TRUE)
proj <- addClusters(input = proj, reducedDims = "IterativeLSI" , force = TRUE)
proj <- addUMAP(ArchRProj = proj, reducedDims = "IterativeLSI" , force = TRUE)
proj@cellColData$Clusters <- ifelse( proj@cellColData$Clusters == "C1" , "C1" , "C2" )

p2 <- plotEmbedding(ArchRProj = proj, colorBy = "cellColData", name = "Clusters", embedding = "UMAP")

out_file <- paste0(out_path, "/plotEmbedding.pdf")
pdf(out_file)
print(p2)
dev.off()

cell_file <- data.frame( cell = rownames(proj@cellColData) , cluster = proj@cellColData$Clusters)
out_file <- paste0(out_path, "/sperm_atac_cluster.tsv")
write.table( cell_file , out_file , row.names = F , sep = "\t" , quote = F )

###########################################################################################
## 计算motif活性的差异
motif_mat_num <- b
#colnames(motif_mat_num) <- sapply( strsplit(colnames(motif_mat_num) , "_" , 1) , "[" , 1  )
motif_mat_num <- motif_mat_num[,rownames(proj@cellColData)]

## 每个细胞分布算peak和motif的激活程度
sco_atac <- sapply(unique(unique(proj@cellColData$Clusters)),function(x){
    print(x)
    sapply(unique(rownames(motif_mat_num)),function(y){
        mean(
            as.numeric(
                as.vector(motif_mat_num[y,which(proj@cellColData$Clusters==x)])
            )
        )
    })
})

sco_exp2 <- data.frame(sco_atac)
sco_exp2$Motif <- rownames(sco_exp2)
sco_exp2 <- sco_exp2[,c( ncol(sco_exp2) , (1:ncol(sco_exp2)-1) )]

out_file <- paste0( out_path , "/MotifScore.sperm.MeanByCluster.tsv" )
write.table( sco_exp2 , out_file , row.names = F , sep = "\t" )

###########################################################################################
## 分别鉴定peak-gene
corrCutoff <- 0.5       # Default in plotPeak2GeneHeatmap is 0.45

p2gMat <- plotPeak2GeneHeatmap(
  proj, 
  corCutOff = corrCutoff, 
  groupBy="Clusters",
  nPlot = 1000000, returnMatrices=TRUE, 
  k=25, seed=1)

kclust_df_out <- cbind( kclust=p2gMat$ATAC$kmeansId , p2gMat$Peak2GeneLinks )
tmp_peak_dat <- data.frame(atac_proj@peakSet)
tmp_peak_dat$peak <- paste0( tmp_peak_dat$seqnames , ":" , tmp_peak_dat$start , "-" , tmp_peak_dat$end )
kclust_df_out <- merge( kclust_df_out , tmp_peak_dat , by = "peak" )

out_file <- paste0(out_path, sprintf("/peakToGeneHeatmap_LabelClust_k%s.tsv", nclust))
write.table(kclust_df_out , out_file , row.names = F , sep = "\t")
